# =============================================================================
# APPENDIX O TABLE 15A: FALSE POSITIVE ANALYSIS (KIM & SUDDUTH)
# Analysis of symbolic legislatures (false positives) effects
# =============================================================================

# --- PACKAGES ----------------------------------------------------------------
library(dplyr)
library(haven)
library(sandwich)
library(ggplot2)
library(reshape2)
library(lmtest)
library(stargazer)

# --- DATA LOADING ------------------------------------------------------------
data <- read_dta("kim_sudduth.dta")
est <- read.csv("estimates_independent.csv")

# --- DATA PREPARATION --------------------------------------------------------
# Ensure consistent identifiers
data$COWcode <- data$cowcode

# Merge and create lagged variables
data_new <- inner_join(est, data, by = c("COWcode", "year")) %>%
  arrange(COWcode, year) %>%
  group_by(COWcode) %>%
  mutate(
    dyn.estimates_lag = lag(dyn.estimates, n = 1),
    dyn.up_lag = lag(dyn.up, n = 1),
    dyn.lo_lag = lag(dyn.lo, n = 1)
  ) %>%
  ungroup()

# Create exact sample for analysis
data_new_exact <- data_new %>% filter(!is.na(dyn.estimates_lag))

data_used_in_models <- data_new_exact %>%
  filter(
    complete.cases(
      dyn.estimates_lag,
      anycoupatt_fin,
      gwf_leadermil,
      ll_rgdpe_PW,
      lag_growth_rgdpe_PW,
      ll_pop_PW,
      postcold,
      anyattyrs
    )
  )

# Create false positive indicator
data_used_in_models <- data_used_in_models %>%
  mutate(
    dyn_quartile_lag = ntile(dyn.estimates_lag, 4),
    false_positive = ifelse(gwf_legislature == 1 & dyn_quartile_lag == 1, 1, 0)
  )


# --- HELPER FUNCTIONS --------------------------------------------------------
cluster_se_vector <- function(model, cluster_vec) {
  sqrt(diag(sandwich::vcovHC(model, type = "HC1", cluster = cluster_vec)))
}

# --- ANALYSIS ----------------------------------------------------------------
# New models with latent measure
m2new <- glm(
  anycoupatt_fin ~ false_positive + gwf_leadermil + ll_rgdpe_PW +
    lag_growth_rgdpe_PW + ll_pop_PW + postcold + poly(anyattyrs, 3),
  data = data_used_in_models,
  family = "binomial"
)

m5new <- glm(
  reshuffleatt_fin ~ false_positive + gwf_leadermil + ll_rgdpe_PW +
    lag_growth_rgdpe_PW + ll_pop_PW + postcold + poly(shuffattyrs, 3),
  data = data_used_in_models,
  family = "binomial"
)

m8new <- glm(
  rechangeatt_fin ~ false_positive + gwf_leadermil + ll_rgdpe_PW +
    lag_growth_rgdpe_PW + ll_pop_PW + postcold + poly(rechattyrs, 3),
  data = data_used_in_models,
  family = "binomial"
)


# Calculate clustered standard errors
cl_ks_merged <- data_used_in_models$COWcode

se_m2n <- cluster_se_vector(m2new, cl_ks_merged)
se_m5n <- cluster_se_vector(m5new, cl_ks_merged)
se_m8n <- cluster_se_vector(m8new, cl_ks_merged)

# --- OUTPUT ------------------------------------------------------------------
stargazer::stargazer(
  m2new, m5new, m8new,
  type = "latex",
  se = list(se_m2n, se_m5n, se_m8n),
  star.cutoffs = c(0.1, 0.05, 0.01),
  notes.append = FALSE
)
